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pi I We investigate the numerical performance of the regularized deconvolution closure introduced recentiy 

by the authors. The purpose of the closure is to furnish constitutive equations for Irwing-Kirkwood-NoU 
procedure, a well known method for deriving continuum balance equations from the Newton's equations 
of particle dynamics. A version of this procedure used in the paper relies on spatial averaging developed 
^ by Hardy, and independently by Murdoch and Bedeaux. The constitutive equations for the stress are 

—i given as a sum of several operator terms acting on the mesoscale average density and velocity. Each term 

is a "convolution sandwich" containing the deconvolution operator, a composition or a product opera- 
tor, and the convolution (averaging) operator Deconvolution is constructed using filtered regularization 
methods from the theory of ill-posed problems. The purpose of regularization is to ensure numerical 
stability. The particular technique used for numerical experiments is truncated singular value decomposi- 
tion (SVD). The accuracy of the constitutive equations depends on several parameters: the choice of the 
averaging window function, the value of the mesoscale resolution parameter, scale separation, the level 
• of truncation of singular values, and the level of spectral filtering of the averages. We conduct numerical 

experiments to determine the effect of each parameter on the accuracy and efficiency of the method. Par- 
tial error estimates are also obtained. 
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1. Introduction 

Many biological and advanced man-made materials are highly heterogeneous, with internal structure 
spanning a broad range of spatial scales. Micron-to-nanometer scales seem to be of particular impor- 
tance for understanding the response of such materials. Traditional continuum models miss important 
effects at sub-micron scales: micro-instabilities, dispersive energy transport, radiative damping, and 
presence of phonon gaps (see e.g. Charlotte & Truskinovsky (2012)). To capture these features, one 
could use molecular dynamics (MD), but analysis of such models is difficult, and their computational 
efficiency is poor. For millimeter size samples, MD simulations can only access time scales on the order 
of 10^** sec, which is unsatisfactory for many applications. This makes it useful to look for mesoscopic 
models having the efficiency of classical continuum and capable of describing smaller-scale effects. 
A related problem in computational science is the design of fast multiscale algorithms for simulating 
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coarse scale particle dynamics. Another relevant question, periiaps the most fundamental of the three, is 
to find a general method for generating continuum constitutive equations from the underlying atomistic 
models. 

Continuum balance equations of can be derived by averaging particle dynamics. This approach was 
originated by Irving & Kirkwood (1950), who used ensemble averages to derive hydrodynamics equa- 
tions directly from Hamiltonian dynamics. Soon after, Noll (1955) improved the mathematical foun- 
dation of the method by introducing finite size "window functions" instead of delta functions. Several 
decades later. Hardy (1982) and Murdoch & Bedeaux (1994), Murdoch & Bedeaux (1996), Murdoch 
& Bedeaux (1997), Murdoch (2007) independently developed a theory based entirely on space-time 
averaging. The averages can be defined for all spatial and temporal scales, and at each scale, no matter 
how small, they satisfy exact continuum balance equations. In that sense, the theory works the same 
way at all scales. This is a useful property because it ensures seamless integration of different regions in 
coupled multiscale simulations. Recently, Hardy-Murdoch approach has attracted much attention in the 
applied mathematics, materials science, and physics communities. It has found use both for solids by 
Tadmor & Miller (201 1) and fluids far from equilibrium by Evans & Morriss (2008). These books also 
contain a large number of references. A few representative examples are E et al. (2009) (complex fluids), 
Lehoucq & Sears (201 1) (peridynamics), and Zinomerman et al. (2010) (solids). In the latter work, the 
authors develop a material (Lagrangian) frame approach as opposed to spatial (Eulerian) formulation of 
the standard Hardy-Murdoch equations. The recent papers by Admal & Tadmor (2010, 2011) extend 
the averaging approach to multibody potentials. 

Despite its many attractive features, Hardy-Murdoch averaging does not produce an effective con- 
tinuum theory. The fluxes in their balance equations are functions of particle positions and velocities. 
This means that trajectories of all particles must be known before fluxes can be calculated. In contrast, 
constitutive equations of classical continuum express fluxes in terms of density, velocity, deformation, 
and temperature. To recover the efficiency of continuum description, one needs to formulate constitu- 
tive equations for Hardy-Murdoch fluxes, that is, approximate them by operators acting on the averages. 
These approximations would play the same role as constitutive equations of classical continuum: they 
would decouple balance equations from the underlying particle dynamics. 

Meso-scale continuum theories differ from classical continuum. The two features they exhibit most 
often are non-locality and scale dependence. For example, spatially non-local solid mechanics theories 
were proposed by Eringen (1976) and Kunin (1982). Charlotte & Truskinovsky (2012) introduced a 
temporally non-local (history dependent) model of one-dimensional lattices with linear nearest neighbor 
interactions. Another non-local and scale dependent approach is peridynamics proposed by Silling 
(2000). In the past 12 years peridynamics has gained popularity and has been the subject of many works 
in applied mathematics, notably by Du, Gunzburger, Lehoucq, Silling and many others (see the review 
papers by Silling & Lehoucq (2010); Lehoucq & Sears (201 1); Du et al. (2012) and references therein). 
The original theory is phenomenological, so the constitutive equations are postulated or fitted from 
experiments. Scale dependence is introduced by using a fixed size window function (called horizon by 
Silling). The latter feature is similar to Hardy-Murdoch procedure, so one expects that peridynamic 
equations can be derived by averaging atomistic models. This was recently done by Lehoucq & Sears 
(2011). They provided an exact description of the peridynamical force state in terms of the particle 
positions and velocities. 

A method for generating non-local and scale-dependent constitutive equations for Hardy-Murdoch 
averages was introduced in Panchenko et al. (2011) and tested numerically in Panchenko et al. (subm). 
Applications to discrete models of fluid flow were studied in Tartakovsky et al. (2011), Panchenko & 
Tartakovsky (subm). Recently, similar methods have been successfully used in large eddy simulation 
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of turbulence by Berselli et al. (2006); Kim et al. (2012); Layton & Rebholz (2012). An outline of the 
method is as follows. The fine scale model is a system of Newton's ODEs for particles. A typical 
interparticle distance is characterized by a small parameter e = N^^l^, where d is the dimension of the 
physical space. The spatial averages are defined using a window function depending on the mesoscale 
resolution parameter 77 e. 

(i) Rewrite non-linear averages as (linear) convolutions of the window function and certain 
functions of particle trajectories. The number of these micro-scale dynamical functions 
is small and equals the number of continuum quantities in the theory. A proper choice of 
the window function insures invertibility of the convolution. Thus the above micro-scale 
dynamical functions are recoverable, at least in principle. 

(ii) The difficulty here is that deconvolution problem is unstable (ill-posed), so small errors 
in the averages may produce large errors in the recovered functions. Instability can be 
overcome by using regularization methods. Regularization produces stable approximations 
of the micro-scale quantities in terms of the averages. 

(iii) In the simplest case of nearly isothermal dynamics, the energy balance is assumed 
to be trivial. The only flux that requires closure is the stress in the momentum balance 
equation. The observable averages in this case are density and linear momentum 
p^v^ . The corresponding recoverable functions are, respectively, Jacobian J of the inverse 
fine scale deformation map, and the product Jv of this Jacobian and a fine scale velocity 
interpolant v. More details on the definition of these quantities can be found in Panchenko 
et al. (201 1, subm) and in the brief description provided in Section 2. 

(iv) Next, we take advantage of a remarkable property of Hardy-Murdoch averages: the 
exact stress can be shown to depend only on the recoverable functions J and v. Thus we 
can substitute their deconvolution approximations into the exact flux equations and obtain 
(approximate) constitutive equations. This produces a meso-scale system in closed form 
that is completely decoupled from the underlying MD. 

The constitutive equations express stress in terms of the average density and average velocity. Fine 
scale information is also incorporated, but the nature of this information is such that solving the ODEs 
is no longer required. Typically we would use the initial conditions and the equation for the interatomic 
potential, but nothing else. The main structural unit of the stress equations is a "convolution sandwich" 

R^SQ^, (1.1) 

where Rr^ is the convolution (averaging) operator, 2^ is the deconvolution operator (this is an approx- 
imate regularized inverse of Rr{), and 5 is a nonlinear composition or a product type operator. The 
complete stress is a sum of the terms of the form (1.1) where \p^\ and Qjj [p''v''] may be present. 

In this paper we study algorithmic realization of the regularized deconvolution closure. Since the 
problem at hand is nonlinear, the error estimates that we provide are not likely to be tight, so we also 
study the error evolution numerically. The objective is to understand the dependence of the error on 
various parameters of the method. Specifically, we study the effects of a choice of a window function i/a, 
resolution parameter 7], scale separation, truncation level in SVD and filtering of spectral coefficients. 
The numerical experiments are conducted as follows. We solve the system of particle dynamics ODEs 
and compute all relevant averages. Next, we apply regularized deconvolution to computed average 
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density and momentum to obtain approximations of / and v. These approximations are used to calculate 
the approximate stress that is then compared to the directly computed ("exact") stress. For simulation 
of particle dynamics, we use two sets of the test initial conditions. The initial positions for both of 
them are uniformly spaced. The initial velocities in the first test case are a discretized low frequency 
sine function. In the second case, the initial velocities are prescribed by a truncated 4th polynomial 
imitating a Gaussian. The corresponding Fourier spectnun is full (unlike the first test case) and decays 
at a polynomial rate. 

A discretization of the integral convolution equation is a linear system. As a consequence of ill- 
posedness, the matrix of this system is ill-conditioned. To regularize the discrete problem, we use 
truncated singular value decomposition (SVD): the exact solution is approximated by its projection 
onto the subspace spanned by singular vectors corresponding to larger singular values. The number of 
the retained singular values plays the role of regularization parameter Since the convolution kernel is 
dynamics-independent, SVD should be computed only once. This can be done prior to running time- 
dependent simulations. The same kernel may be used with different dynamical systems, and thus the 
cost of computing SVD can be excluded from the operation count of the dynamical simulation. The 
algorithms for computing SVD are readily available from the standard Unear algebra packages such as 
LAPACK Anderson et al. (1999). All of the above makes SVD-based regularizations convenient and 
easy to implement. However, using SVD has several drawbacks. Many standard SVD algorithms are 
iterative, and smaller singular are computed less accurately. The same is true for computation of the 
associated singular vectors. In addition, the performance of the standard algorithms including LAPACK 
becomes worse as the condition number of the matrix increases. For numerical experiments performed 
in this work, advantages of using the truncated SVD seem to outweigh the disadvantages. 

In the continuum case, degree of ill-posedness depends on the rate of decay of the singular values of 
the kernel (see Kirsch (1996)). This rate in turn depends on the smoothness of the kernel (see Hansen 
(1987)). It is then reasonable to expect that the rate of decay of singular values plays a similar role 
in the discrete setting. However, the overall error, that is, the error in computing the approximate 
stress (1.1), does follow this pattern. The Gaussian kernel, which has the fastest decay rate among the 
considered window functions, corresponds to a smaller overall error than less smooth piecewise linear 
kernels. We consider six different window functions: characteristic function of an interval, piecewise 
linear trapezoid function and piecewise linear triangular function (finite element function), truncated 2nd 
and 4th order polynomials and Gaussian function. Among these functions, the Gaussian and truncated 
4th order polynomials give the least error. Occasionally the tnmcated 4th order polynomial produces 
slightly smaller error than the Gaussian but using the Gaussian requires a significantly smaller number 
of singular values to achieve comparable accuracy. 

An important aspect of the method's performance is dependence of the error on scale separation 
characterized by, for example, the product of tj and A'^. Ideally, the error should decrease as scale 
separation increases. We test the method in two regimes: one corresponds to increasing 77 keeping N 
fixed, and the other consists of increasing A'^ while tj is fixed. These regimes are not equivalent. The 
first corresponds to taking a fixed fine scale dynamical system and increasing the size of the averaging 
window. The second regime is different: we keep the size of the window fixed, but the underlying 
particle systems vary as follows. We increase A'^ while keeping the total mass fixed, and rescale the forces 
so that the total energy E (N) is bounded independent of A'^. The initial conditions in this regime are given 
by increasingly finer discretizations of the same continuum functions. Existence of the uniform in A' 
bound on energy does not imply monotonic dependence on A'^, and thus it is unlikely that the averages 
of different particle systems depend on A^ monotonically. 

In the first regime, as tj increases, the density and Unear momentum tend to get smoother and 
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converge to zero as 7] oo (see equations (6.3) and (6.4)). The exact stress depends on 77 similarly, 
but the dependence of the approximate stress on 77 is more complicated. To form the approximate 
stress, we replace the exact recoverable functions by their deconvoluted approximations. The latter 
depend on rj, and thus the nonlinear functions of the recoverable functions appearing in (1.1) become 
7] -dependent. Finally, to produce the approximate stress, these nonlinear functions are averaged with the 
TJ -dependent window function. The quality of the approximation is thus determined by two competing 
tendencies. On the one hand, for larger 77, one expects worsening of the reconstruction Qrj because the 
same amount of fine scale information is harder to extract from increasingly smoother and smaller data. 
On the other hand, the averaging in (1.1) decreases the total error by damping the high frequency 
modes. The error in the deconvolution is primarily at high frequencies. Application of the nonlinear 
dispersive operator S may cause propagation of the error into low frequency harmonics. This effect 
may be difficult to counteract by applying but in simulations the overall relative error tends to 
decrease with increasing 77, even though the deconvolution error increases with 77. This seems to be 
a remarkable feature of the deconvolution closure: the stress is computed more accurately despite the 
fact that the recoverable functions are approximated more poorly. This shows that the constitutive 
equations introduced in Panchenko et al. (201 1, subm) become more accurate at larger scales. We also 
confirm this by comparing Fourier spectra of the exact recoverable functions and stresses to the spectra 
of their approximations. For both test cases, the truncated SVD approximation captures low frequency 
harmonics but looses high frequencies. The Fourier spectra of the stress are dominated by low frequency 
harmonics and they are approximated well. It seems that dispersion due to the presence of 5 in (1.1) is 
relatively weak. 

In the second regime, we simulate evolution the error as N increases. We observe that for certain 
times the error is not monotonic in N. This deviation from monotonicity seems to be related to the 
behavior of the computed total energy of the system. While the exact total energy is conserved, the 
computed energy oscillates in time in a nearly independent of manner. The oscillations are between 
the initial value (the exact energy of the system), and a shghtly larger (by at most 0.05%) value. At 
times when the computed energy is close to the initial value, the error decreases as A'^ increases. At 
other times, the dependence of the error on is not monotonic. It appears that enforcing conservation 
of energy numerically, for example, by using energy-preserving Runge-Kutta methods (see Celledoni 
et al. (2009)), should ensure monotonic decrease of the error with A^. 

The classical theory of ill-posed problems Kirsch (1996); Tikhonov & Arsenin (1987); Morozov 
(1984) deals mostly with operators between Hilbert spaces, because for integral operators, the Hilbert 
space structure is required for existence of singular value decomposition. In the discrete setting, one 
can use SVD with any p-nona, including the 00-norm. In the paper, we provide such error estimates for 
p G [1,°°]. The values of p for the right hand side and the solution may be different. We also use the 
00-norm for computational assessment of the quality of approximation. 

The organization of the paper is as follows. The background concerning the fine scale problem 
and spatial averaging is briefly described in Section 2. Section 3 contains the necessary facts about 
filtered regularization methods in the discrete setting. In Section 4 we study various window functions 
used to set-up averaging and the effect of the choice of a window function on the quality of the stress 
approximation. In Section 5 we analyze the effect of the resolution parameter on the reconstruction of 
microscopic quantities of interest and eventual stress approximation. In Section 6 we study the Fourier 
spectra of the reconstruction error and the overall error in the stress. The effect of the scale separation is 
considered in Section 7. Error estimates for filtered regularization methods for /? G [1 , 00] are provided in 
Section 8. Finally, in Section 9 we give conclusions. In appendices, we define various window functions 
used in the numerical experiments and provide the formula for the Lennard- Jones potential. 
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2. Governing equations 

2. 1 Micwscale equations 

The object of our study is a system of particles that move according to Newton's law of motion, 
where N^l. Particle positions and velocities are denoted by and V;(f), i= 1 , . . . ,A^, respectively. 
For simplicity, we assume that all particles have the same mass M/N where M is the total mass of the 
system. The equations of motion are 

Qi = Vi, (2.1) 
^v, = + (2.2) 

with initial conditions 

9,(0) =9?, v,(0)=vO. (2.3) 

Here ff^'^ is an external force (gravity) and /, = Y,j fij^ where fij is an interaction force exerted on 
a particle i by a particle j. These forces are generated by a pair potential U{\qi — qjl). In numerical 
experiments we use the classical Lennard- Jones potential (see equation (A.l)). 

To study the behavior of the system for large A^, it is convenient to introduce a small parameter e = 
N^^/'^, where c/ = 1 , 2, 3 is the dimension of the physical space, and a mesoscopic resolution parameter 
T] satisfying < t] < 1 and e <C i]. Using these parameters we define the microscopic length scale 
eL and the mesoscopic length scale tjL, where L = diam(i2) and Q is the computational domain. To 
make the total energy of the system bounded independent of N, we scale interparticle forces by e as in 
Panchenko et al. (subm). 

2.2 Mesoscale averaging and dynamics 

The first step of model reduction in Panchenko et al. (2011, subm) is to use space-time averaging pio- 
neered by Noll (1955), Hardy (1982) and Murdoch & Bedeaux (1994); Murdoch (2007). In this section, 
we briefly describe the averaging approach to make the exposition self-contained. For simplicity we con- 
sider only spatial averaging and refer to Noll (1955), Hardy (1982) and Murdoch & Bedeaux (1994); 
Murdoch (2007) for more details. 

To set up spatial averaging, choose a fast decreasing window function ^f{x)^0 such that /TL W{^) = 
1. We assume that \^ is continuous and differentiable almost everywhere on the interior of its support. 

After scaUng by scaling by tj , the function ^/^((x) = ^i/A is used to define the averages of micro- 
scopic quantities. For example, the approximate number of particles within the distance tjL of a point x 
is 

N 

n'^{t,x) = Y,Vri{x-qi{t)). 
Similarly, one defines the average density and linear momentum: 

P''(^^) = f fv^nl^-^iW), (2.4) 

i=i 

r^^(t,x) = ^l,Vi{t)Wn{x-qM. (2.5) 
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Differentiating equations (2.4), (2.5) with respect to t, and using the ODEs (2.1), (2.2) yields exact 
balance equations for and p^v^ (see Hardy (1982); Murdoch & Bedeaux (1994)): 

^,p''+div(p''v'')=0, (2.6) 

a,(p''v'')+div(p''v''(g)v'')-div7"' =0. (2.7) 

Here the system is assumed to be isolated (/-^^ =Q)\T^ is the total stress written as = Tj^^ + T^^^^ 
(see Murdoch (2007)), where 

T'l'^) {t,x) = - £m,(v,- - v'' (f )) (v; - v'' V^^ {x - q^) (2.8) 
is the convective stress, and 

T'^{t,x)i^i,,) = Y.fij®{qj-qi) Cwn {s{x-qj) + {\-s){x-q,))ds (2.9) 

(<■,;) •'^ 

is the interaction stress. The summation in (2.9) is over all pairs of particles 7) that interact with each 
other. A similar energy balance equation can also be written if needed (see Hardy (1982); Murdoch 
& Bedeaux (1994)). We refer to (2.6), (2.7) as the meso system. The unknowns p^ and are at the 
mesoscale whereas the stress 7"' depends on microscopic variables ^, and v,. Hence, the meso system 
is not closed because one must know the trajectories of all particles before can be determined. 
Therefore, solving the exact meso system is no more efficient than solving the original ODEs (2.1), 
(2.3). 



2.3 Integral approximations and closure 

To improve efficiency, we use the closure method introduced in Panchenko et al. (201 1, subm). In this 
section, we outline the method for convenience of the reader. 

The objective is to find an approximation of in terms of p^ and p^v^. First, observe that the 
discrete sums in (2.4), (2.5) can be approximated by the convolution integrals 

Rr, [/] (x) = Jvri{x- y)f{y)dy « . (2.10) 

Here Rr^ is a convolution operator, / is a microscopic quantity to be recovered, / is the known 
mesoscale average. The integral approximations of the average density and average momentum are 
as follows. 

M Ji Mr 
= ^Ly'nix~9i{t))^j^J^Wri{x-q{t,X))dX (2.11) 

M f M 

M Mr M 

p'^v^{t,x) = -I2v,(r)r^(^-9,(f)) -J^J^ Wn{x-y)v{t,y)J{t,y)dy = — 7?^[v7](x). (2.12) 
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Here M is the mass of the system, and |^2 1 is the volume (Lebesgue measure) of the domain ^2 occupied 
by the system. The functions q{t,X), v{t,q) are suitable interpolants of particle positions and velocities. 
The Jacobian 

/=|detV4-'| (2.13) 

describes local volume changes. 

Integral approximations of the stress are obtained similarly (see Panchenko et al. (subm)). To 
approximate the interaction stress T^-^^^ we write = q(t,X\),qj = q{t,X2), and approximate the 
double sum in (2.9) by the double integral with respect to Xi , 2. After several changes of variables of 
integration we arrive at 

Equation (2.14) can be written compactly as 

where Rri is the convolution operator, and the operator Si is a composition of shifted multiplication 
J{t,R) J{t,R+ ^p)J{t,R- ^p) and p-integration with the weight \£2\~'^\p\~^p iS> p U'{\p\). 
Convective stress T^^^ can be approximated similarly: 

Tl){t,x) « -^J^Wvix-ymt,y)-v'^{t,x,))^{vit,y)-v'^{x,t))J{t,y)dy (2.16) 

= RnS2[J,v], (2.17) 

where the operator ^2 is a multiplication of — by the diadic product of v{t,y) — {t,x) with itself. 

Equations (2.14 and (2.16) show that stress can be represented as an operator acting on the two fine 
scale functions J and v. These functions can be approximately recovered by inverting the convolution 
operator Rrj (see (2.11) and (2.12)). The difficulty here is that straightforward inversion is unstable. 
Indeed, 7?,) is compact, hence it is not surjective. This implies that ' is unbounded, and thus small 
perturbations of the right hand side may produce large perturbation in the computed solution. On the 
discrete level, the matrix discretization of Rri will have a large condition number. Therefore, the decon- 
volution problem (reconstructing / from the knowledge of R^j [/]) is ill-posed. Such problems can be 
approximately solved using regularization (see, e. g. Kirsch (1996); Morozov (1984); Tikhonov & 
Arsenin (1987)). Regularization consists in approximating the exact unstable problem by a parameter- 
dependent stable problem. The solution operator Qrj of this problem provides an approximation to 
the exact inverse operator. Computational implementation of regularized deconvolution is described in 
Section 3 below. 

The meso system (2.6), (2.7) is closed by inserting the approximations 

J^M-^\Q\Qr,n v-^^ (2.18) 

into (2.15), (2.17). The resulting operator structure is a sum of terms of the form (1.1), acting on and 
p''v''. 

Computational realization of the closure equations is obtained by approximating microscopic posi- 
tions using (2.13). In one dimension, this is particularly simple, since the Jacobian coincides with the 
gradient. Approximating it by standard finite difference discretizations, we can generate approximate 
positions. Then approximate positions and velocities are inserted in (2.8) and (2.9). 
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2.4 Test cases 

For computational testing we use Lennard-Jones potential defined in (A.l) and two sets of the initial 
conditions. The problem is assumed to be periodic with period L. The initial positions in both cases are 
equally spaced with 



q^j ^ yj - -J Ax, where ; = A?. (2.19) 

First test case. The initial velocity is a one mode sine function 

v(''(ji:,0) = 10^2 sin x L. (2.20) 

Second test case. The initial velocity is a continuous function that is a truncated 4th order polynomial 
on [| , ^] with double roots at x = j and x = | and zero otherwise: 

y(2)(x,0) = <^ ^ ^ ^ ^ (2.21) 

I 0, otherwise. 

The initial velocities are shown in the left panel of Fig. 1. The main difference between these initial 



1\ . . . L 
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0.01 
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Fig. 1. Initial velocities and their Fourier spectrum. Left panel shows initial velocities. Red solid curve is for the first test initial 
velocity v'''(a',0) defined in (2.20), black curve - second test initial velocity i''-'(a',0) given in (2.21). Right panel has discrete 
Fourier transform of v'-' (x, 0) . 



velocities can be seen by looking at their Fourier spectra. The first initial velocity has only one low 
frequency Fourier mode while the second initial velocity has all Fourier modes, shown in the right 
panel of Fig. 1, present that level off at 10^^ for = 1000. As increases, the tail of the Fourier 
spectrum levels off at a slightly lower value (not shown here). For example, for = 10,000, the tail is at 
10^'^. Therefore, we expect to see more complicated nonlinear dynamics and a stronger effect of high 
frequencies in the second case. 

The ID version of system of ODEs (2.1), (2.2), (2.3) with the initial conditions defined in (2.19), 
(2.20) and (2.21) is solved using the Velocity Verlet method until t — \ for various A^. 

3. Filtered regularization methods 

On a discrete level, equations (2. 1 1), (2. 12) reduce to a linear system 



(3.1) 
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where ^7 is a known average quantity such as -jfP^ and -jf'p^v^ and x is either J or vJ. To produce 
b, averages are discretized on a coarse mesh with B nodes, B ~ 1 /t] and B <^N, and the solutions are 
rendered on a finer mesh with N' nodes. In practice, A'^' can vary between B and A'^. In our numerical 
experiments we use A^' = A^. 

The matrix , obtained by discretizing the kernel V^,j , has dimensions BxN' and r = rank(A'' ) ^B. 
There are two difficulties associated with solving (3.1): (i) the condition number of is large because 
of the ill-posedness; and (ii) the system is underdetermined and has multiple solutions. 

Performing singular value decomposition (S VD) of A^ we obtain r non-zero singular values Oj and 
left and right singular vectors and ^j, satisfying 

A^j^Gj^j, {A^f^j = aj^j, 7 = l,2,...,r. (3.2) 

The singular vectors | j, |^ have length B and A^, respectively, and are orthonormal in their respective 
spaces. In this paper we work only with matrices A^ satisfying 

(7,e(0,l]. (3.3) 

This assumption holds for discretizations of all convolution kernels under consideration. 
Expanding the right hand side 

we can write the minimal norm solution of (3.1) 

^=Y,^j^P ^J = zr' y = l>2,...,r. (3.5) 
j=i 

This solution is orthogonal to the null-space of A^ . When the condition number of A^ is large, the 
solution is highly unstable. To stabilize the computation, one can use regularization (see, e. g. Kirsch 
(1996); Hansen (1987)). Here we limit ourselves to filtered linear regularization methods that replace 
the exact solution (3.5) with the approximation 

.« = I^,^^|,. (3.6) 

j=i 

The function is called a filter function. Generally it should satisfy (see Kirsch (1996)) 

1. |<^(i,a)|<lforalla>Oandie(0,l]; 

2. for every a > there is c(a) such that |(^) (s, a) | ^ c{a)s for all s G (0, 1] ; 

3. \\ma^o<^{s.a) = 1 for every s e (0, 1]. 

By Theorem 2.6 in Kirsch (1996), (3.6) defines a regularized approximate inverse Ra to A^A that has 
the property ^raa^oRa{A^A)x^x'm 2-norm for each jc e R^. 
Examples of filtered methods are Tikhonov regularization with 
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truncated SVD with 



a 



{ 



1 if Uj < Oa 
otherwise, 



(3.8) 



and Landweber iteration with 



.2\n+\ 



(3.9) 



In (3.8), the role for regularization parameter is played by the cut-off value Oa, while for Landweber 
iteration the regularization parameter is the reciprocal of the number n of iterations. The SVD approach 
may be more efficient than iterative methods (see Hansen (1987)) because SVD can be pre-computed 
and used at each time step for different right hand sides. The algorithmic realizations of (3.7)-(3.9) and 
other numerical regularizing techniques are discussed in the book by Hansen (1987). 

In this paper, we use a truncated SVD method to compute the minimum norm solution (3.5) and 
we retain only those singular values whose magnitude is greater than = 10"'^. In addition we 
also use spectral filtering technique on the right hand side b, which is similar to the Fourier filtering 
used in Krasny (1986). The coefficients bj that lie below the tolerance value tol = 10^'^ are set to 
0. We find this filtering helpful when both Uj and bj are close to the machine precision. In that case, 
both bj, and xj may carry a large relative error. Better reconstructions are obtained by discarding the 
corresponding coefficients x,- = bj/Oj in the computed solution. Varying tol between 10^'^ and 10^'^ 
does not essentially change approximations to stresses, while setting the filter level at a higher value 
increases the error as expected. 

To compute SVD of A'', we use a standard LAPACK subroutine DGESVD. This routine computes 
singular values and left and right singular vectors. It should be noted that the accuracy of this subroutine 
decreases as the singular values decrease. Specifically, each computed singular value a, differs from 
true value (7, by at most 



where p{m,n) is a modestly growing function of m and n, e = 2^ « 1.1 1 • 10^ is machine epsilon, 
and (7i is the largest eigenvalue. Thus singular values near Oi are computed to high relative accuracy and 
small ones may not be. The computed singular vectors are always orthogonal to working precision but 
their accuracy depends on how close they are to each other, i.e. if (7, is close to nearby singular values. 
In the present case, small singular values are densely spaced, so the corresponding singular vectors may 
be inaccurate (see Anderson et al. (1999)). 

Jacobi's method (see Demmel et al. (1999); Demmel (1997); Demmel & Veselic (1992); Slapnicar 
( 1 992)) is another algorithm for computing SVD. Unlike conventional algorithms, this method is capable 
of high relative accuracy even for smallest singular values. The original version was slower than the 
standard algorithms but the speed was recently improved by Drmac & Vesehc (2007a,b). Subroutines 
DGEJSV and DGES VJ are available as a part of the current version of LAPACK. It would be interesting 
to see if using Jacobi's method improves the accuracy of the stress approximation. 

4. Choice of a window function 

We restrict our attention to the window fimctions ^ satisfying the conditions 

\j/ is non-negative, continuous, and differentiable almost everywhere on the interior of its suppor(4dd 



(7,| ^ p{m,n) • e • (7i 




(4.2) 



Va(x) ^0 as |x| oo. 



(4.3) 
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The function y/{x) can be compactly supported or fast decreasing as Gaussian. We consider several win- 
dow functions: piecewise constant (characteristic) \j/^^\x), piecewise linear trapezoidal shape \j/^-^\x), 
piecewise linear triangular (x), truncated quadratic \jf^^^ (x), truncated 4th order polynomial i/a*^^^ (x) 
and Gaussian i//f^^(x). The functions are chosen to have the compact support on [—L/2, L/2]. For the 

Gaussian function xj/'^^^ (x) = exp ^— , the standard deviation is (7 = L/6, i.e. only about 0.2% 
of the area under is outside [—L/2, L/2]. Window functions are depicted in Fig. 2. The shifted 









-L/2 L/2 


/ 




\ 




-L/2 -L/6 L/6 L/2 

1 .5/L 1 .875/L 6(2it)-"^/L 





Fig. 2. Window functions. Top left: piecewise constant (characteristic) \l/^^\x); top middle: piecewise linear trapezoidal shape 
(^'^'(jc); top right: triangular function i/(^'(;c); bottom left: truncated quadratic <ff'^*\x); bottom middle: tnmcated 4th order 
polynomial \lf^^\x); bottom right: Gaussian \j/^^\x) 

and rescaled window function 1//^, {x — qi{t)) defines the average properties of the particles in the vicin- 
ity of X. In the case of the characteristic function \ir^^\ the contributions of all particles within a fixed 
distance of x are weighted equally. A better choice of a window function is obtained when the weight 
decreases to zero as the distance between the particle and the observation point increases. This adds 
another "desired" characteristics of a window function (see Root et al. (2003)): 



V/'(x) has a maximum at x = 



(4.4) 



Clearly, functions \if^^\x) and i/aP^ (x) do not satisfy this property, whereas the rest of the functions does. 
A window function can have a different degree of differentiability, ranging from piecewise constant (as 
characteristic function) to infinitely many times differentiable as Gaussian. The order of differentiability 
affects the speed of decay of singular values to zero. In the left panel of Fig. 3 we plot singular values of 
Va(') (x), j = 1 , . . . , 6, for tj = 0.1, at = 1000 and B = 500. Results are similar for N = 10, 000. Since we 
have to invert an operator with kemel , increasing smoothness of the kernel increases ill-posedness of 
the inverse problem. On a discrete level, the situation is somewhat different. First, discretization itself 
regularizes the inverse problem (see Kirsch (1996)). In addition, truncating SVD provides additional 
regularization (see, for example, Hansen (1987)). It is natural to ask how smoothness of the window 
function effects the reconstruction quality. The least smooth function is a characteristic function i/^^' (x), 
followed by piecewise linear trapezoidal y/'^^) (x) and triangular function (/a^^) (x). Then we have a trun- 
cated quadratic i//*^^^(x), truncated 4th order polynomial i/a(^)(x) and Gaussian i/a(^)(x). The singular 
values of the corresponding matrices are given in the left panel of Fig. 3. As expected, the decay of 
spectral coefficients is fastest for the Gaussian and slowest for the characteristic function. A sharp drop 
of singular values for the characteristic function i/a(^)(x) and triangular function \^^^\x) indicates that 
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the discrete problem is numerically rank-deficient. The same can be said about the Gaussian that had 
only about a third of singular values above the machine zero. 

If the convolution kernel y/,] and function / is periodic or extended periodically, then Rrj is a circular 
convolution operator (see Mallat (2009)). If both / and [f] are discretized on the same grid, then the 
eigenvectors of the circular convolution operator are the discrete complex exponentials and the eigen- 
values are Fourier modes of the window function \j/ri (see e.g. Mallat (2009)). Therefore, eigenvectors 
corresponding to the smallest eigenvalues of the Rri carry information about high frequency compo- 
nent of a solution. If R^ [f] and / are sampled at different scales, we have to deal with singular value 
decomposition instead of eigenvalue decomposition. In this case, singular vectors corresponding to the 
largest singular values would have contribution not only from low frequencies, but also from some high 
frequencies. Nevertheless, singular vectors corresponding to smallest singular values will be the most 
oscillatory and we can still think that singular vectors corresponding to the smallest singular values 
represent the oscillatory part of a solution. In this regard, the fact that singular values for all window 
functions but Gaussian decay slowly indicates that solutions corresponding to these window functions 
wiU have high frequency component present. Hence, if there is a numerical error in computation of 
singular values and singular vectors, especially in those corresponding to the smallest singular values, 
there may be a significant error in computation of a solution. On the other hand, singular values of 
a scaled Gaussian window function decay very fast. This means that a fewer singular values can be 
used to represent a solution. Thus, using a Gaussian is more efficient for large systems of ODEs if the 
associated error is comparable with other choices of window functions. 




200 400 600 200 400 600 



Fig. 3. Left panel: singular values for window functions I//''', / = 1, ... ,6, with T] = 0.1, N = 1000 and B = 500. Right panel: 
singular values for for various 77, 0.01 ^ 77 ^ 0.9, !^ = 1000, B = 500. In both cases, the graphs of singular values with 
N = 10000 and B = 500 are similar. 

In order to approximate the stress, one needs to approximate the exact microscopic positions and 
velocities. These approximations are obtained by generating deconvolution approximations of J and 
vJ, recovering microscopic velocities v by dividing vJ by J, and reconstructing microscopic positions 
q from J using (2.13). These are then used to compute approximate stresses. Each of the above steps 
carries some error, and smaller deconvolution error does not necessarily yield to smaller overall error in 
approximating stresses T^^^^^ and r|J,j. 

In Fig. 4 we compare the /"-relative errors in approximation of the convective T^^ and interaction 

'^{int) Stresses for the first test case. Clearly, the characteristic function has the worst performance. 

The absolute error (not shown here) in T^^,^ is at most 10^^ or 13%, while the error in using other 

window functions is much smaller. For example, piecewise linear gives the absolute error of two 
order lower: at most 2.5 • 10^^ or 6%, while the 4th order truncated polynomial \jf^^^ and Gaussian 
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give the error 10 or only 2 — 3%. The error in using the rest of the window functions is between 10 
and 2.5 • 10^^ or between 3 and 6%. The right panel of Fig. 4 suggests that the approximation of Tj^^^^ 

with i/A^'' has too large error, which is not acceptable. The error with y/'^' is at most 4 • 10^^ (reaching 
75 — 100% at some times) and drops to 2 • lO"'* (or 0.5%) with xff^^^ and i/^^'. Functions i/'^' and i/'"*' 
give 1 — 3% error 

The worse behavior of followed by could be explained by the fact that function i//''^ is 
not even continuos everywhere that violates condition (4.1), while i/'^' is only piecewise continuous. 
Moreover, both functions do not have a strict maximum at x = 0, thus violating the condition (4.4). For 
the rest of the window functions considered here, all conditions (4.1)-(4.4) are satisfied. 




Fig. 4. Effect of the choice of the window function on the eiTor in approximation of stresses in the first test case with N = 1000, 
B = 500 and 7] =0.1. Left panel: convective stress. Right panel: interaction stress. 




The results for the second test case are shown in Fig. 5. Again, the worst performance is by xj/^-^K 
followed by y/'^'. The best results are obtained with y/'^' and the rest of the functions give 
intermediate results. Specifically, the absolute error in approximating the convective stress IJ^j using 

V/('' is at most 5-10"'' or 17%, it drops to 1 .5 • 10"^ or 8% with yf^^'> and reaches only 0.5 - 0.75 • 10"*' or 
2 — 5% with y/'^'' — . We also note that the error with functions xjf^^^ — y/'^' is highly oscillatory while 
it is not oscillatory with the rest of the functions. The absolute error in approximation of the interaction 
stress . using function y/'^' is at most 10^^ or 30% (decreases to 15% at later times). Function 
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xffW gives 10^^ absolute error or 4%, while functions \(f^^\ \jf^^'> and \jf^^^ produce only 0.6 — 0.7 • 10^^ 
absolute error or at most 1 — 2%. 

The experiments in this section show that the most accurate approximation of the stresses is obtained 
using either truncated 4th order polynomial function y/'^' or Gaussian The absolute error is typ- 
ically smaller with y/'^' but the relative error is sometimes slightly smaller with Xj/'^^K However, the 
difference in performance of these functions is not significant. At the same time with y/'^' and, for 
example, with = 10,000 and B = 500, one has to use all 500 singular values and singular vectors 
since all singular values are above the threshold Oa — 10^'^ (they decay to at most 10^^, see the left 
panel of Fig. 3), and only 147 singular values with Xj/'^^K which is more efficient for large systems of 
particles. For this reason, in what follows we fix Gaussian window function and study the effect of 
other parameters such as averaging width rj and scale separation. 



5. Choice of mesoscale resolution parameter rj 

The parameter rj in (2.4), (2.5) determines the size of the averaging region and the amount of high 
frequency filtering. For smaller Tj, there is less damping of high frequency content of the solution, 
while using a larger rj produces smoother and smaller averages. The latter is discussed in Section 6. 
As rj increases, the singular values the corresponding matrix decay at a higher rate. This can be 
seen in the right panel of Fig. 3 where we show the singular values for = 1000 and B = 500 and 
rj — 0.01,0.05, ... ,0.9. (The results with = 10000 and B = 500 are similar since they are more 
sensitive to the choice of B and not of A^.) 

On the one hand, for larger T], the average contains less high frequency information, and it may be 
more difficult to reconstruct the same microscopic solution from increasingly smoothed averages. This 
may increase the error in computation of Qri . On the other hand, computation of the stress involves 
another averaging (the outer layer Rj^ in (1.1)) that may decrease the overall error even if J and v are 
recovered more poorly. In this section, we study the cumulative effect of these two competing tendencies 
on the overall relative error in approximating stress. 




0.5 1 0.5 1 



Fig. 6. The effect of the choice of the resolution parameter rj on the error in approximation of the Jacobian (left panel) and 
velocity (right panel) for the first test case with N = 1000 and B = 500. 

To see how the choice of Tj affects the reconstruction of the Jacobian and microscopic velocity, and 
subsequent stress approximation, we fix A^ = 1000, B = 500 and the window function \jf^^^ (x), and vary 
7] between 10"^ and 0.9. 
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5.1 First test case 

In this case, the errors in the Jacobian and velocity approximation, shown in Fig. 6, are essentially 
independent of T], though there is a slight dependence on 77 for times after t = 0.8. Both absolute 
and relative errors in the Jacobian oscillate in a quasiperiodic manner. The amplitude of oscillations 
in absolute error increases slightly with time and reaches 9 • 10^^ at most during the simulation time, 
which is less than 0.01% error. The behavior of the error can be connected to the evolution of the total 
computed energy of the system shown in left panel of Fig. 7. As can be seen from the graph, the energy 
is not completely conserved. Instead it oscillates periodically and deviates from its value at ? = (the 
true energy of the system) by at most 2.5057 • 10~^ (less than 0.05% change). In the regions where the 
energy starts deviating from its initial value, the error in Jacobian reconstruction increases. When the 
energy comes back to its initial value, the error in Jacobian also decreases. The relative error in the 
Jacobian reconstruction is similar to the absolute error since the exact Jacobian has values close to 1. 
The error in velocity reconstruction, shown in the right panel of Fig. 6, has a more nonUnear dynamics. 
Similarly to the error in the Jacobian, the error in velocity increases at those times when the total energy 
starts deviating more from its initial value and decreases when the total energy comes back to this value. 
The maximum absolute error during the simulation time is under 3.5 • 10~^ which is 1.8% error. 



Total energy 



Total energy 



-0.05150 



-0.05152 




-0.05151 



-0.05152 




Rg. 7. Evolution of the total energy for iV = 1000, B = 500 and 17 =0.1. Left panel: the first test case. Right panel: the second 
test case. 



The errors in both convective and interaction stresses depend on 77 as can be seen from Fig. 8. While 
the absolute error in 7"^^^ is not monotonic in Tj at all times, it is typically larger for larger X] . The time 

oscillations of the error do not exceed 1.5 • 10^^ . It is the largest with 77 = 0.8, followed by 77 = 0.6 
and 77 =0.9. However, the relative error, shown in the left panel of Fig. 8, decreases monotonicaUy as 
77 increases. More specifically, the error with 77 = 0.01 is the largest and varies between 9 and 100%, 
whereas it is the smallest with 77 = 0.9 and it varies between 0.3 and 0.6%. For example, for 7j = 0.1, 
the error fluctuates between 0.5 and 3%. Both absolute and relative errors in approximation of T^T^^^ 
oscillate in time and depend monotonicaUy on 77 : they decrease as 77 increases. The latter is shown 
in the right panel of Fig. 8. The absolute error is of the order of 10~^ and reaches 0.7% at most for 
77 = 10~^ and less than 0.5% for the largest 77 = 0.9. 



5 . 2 Second test case 

The simulation results with various 77 are presented in Figs. 9 and 10. Fig. 9 indicates that the recon- 
struction of both Jacobian and velocity gets worse as rj and time t increase. However, the situation 
with stress approximation is different. The error in approximating the interaction stress is the smallest 
with the largest 77 = 0.9 used. It is below 1% during the entire simulation time, whereas the error with 
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the smallest rj — 0.01 reaches 10% by the end of simulations and it is the largest among all values 77 
used. As for the convective stress, the error with rj ~ 0.9 is the smallest at early times (around 0.2%), 
then starts increasing and it is of the same order as with rj = 0.6 and 0.8. After t = 0.6, the errors with 
rj = 0.4 and 0.2 are the smallest (up to 3%) and errors with rj — 0.6, 0.8 and 0.9 reach 4 to 7%. It 
should be noted though that the convective stress is much smaller than the interaction stress (by at least 
two orders of magnitude) and the absolute errors in approximating the convective stress are at the order 
of 10^^ for all Tj used in the experiments. 

Consider, for example, rj =0.1. From Fig. 9 we can see that the error in Jacobian and velocity 
approximation is 0.02% and 2% at early times and reaches 0.12% and 15% by the end of simulations, 
respectively. At the same time, the error in the convective stress varies between 1 — 2% and 5% during 
the simulation period. The error in the interaction stress does not exceed 2%. As we can see, the error 
in the approximation of the Jacobian stays small during the simulation period, though it grows slowly 
with time (at most linearly). Despite of the large error in the approximation of velocity, we still get very 
good stress approximation. The situation is similar for other values of 77 . 




Fig. 9. Effect of the choice of 7) on the error in reconstruction of Jacobian (left panel) and velocity (right panel) in the second test 
case with /V= 1000,8 = 500, 7] =0.01, 0.05, 0.1,..., 0.9. 
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6. Spectral evolution of averages and stresses 

As was mentioned in Section 5, parameter rj determines the size of the averaging window and the 
amount of high frequency fihering. With larger rj, the averages are smoother and smaller This can be 
seen by considering the Fourier transform of an average as follows. 

Recall a typical one-particle dynamical function used in statistical mechanics 

N 

gsm{t,x) = Y,g{q,{t),Vi{t))5{x^qiit)), (6.1) 

;=1 

where 5 is delta-distribution. The Fourier transform of gsm with respect to x is 

gs,nit,^) = l^giq,{t),v,{t))e-^''-. (6.2) 

!=1 

Now compare this with a windowed spatial average 

N 

g{t,x) = Y.8i9,it)Mt))Wnix - qiit)), (6.3) 
and the corresponding Fourier transform 

f(f,x) = f^g{q,{t)Mt)mri^y^'' = Wivm^m, (6.4) 

i=l 

where i// denotes the Fourier transform of Thus, the Fourier transform of g is obtained from gsm by 
low-pass filtering (multiplication by 1/(77^)). If lim|;t|^>>o IV'^('fc)| = 0, which is true for any Li function 
by Riemann - Lebesgue Lemma, then \jf{ri^) converges to for each ^ 7^ as 7] — > 00. Since gsm 
does not depend on rj, equation (6.4) implies g approaches 0. The rate of decay of xjf increases with 
smoothness as xj/. Thus, increasing rj produces progressively more filtered versions of gsm- 

Flexibility afforded by varying rj is convenient for studying large scale behavior of the averages. 
For example, in statistical physics, the derivation of hydrodynamical equations and computation of fluid 
viscosity by Green- Kubo formulas (see e. g. Berne (1977)) employs truncated Taylor expansions of gsm 
at ^ = 0. Estimating the error of these approximations for large ^ may be difficult, while multiplication 
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Fig. 11. Discrete Fourier coefficients of the Jacobian (left panel) and velocity (right panel) for the second test case with N = 
10, 000, B = 500, T] =0.1. The logarithm of coefficients' amplitudes is plotted against wavenumber A at f = 0.9. The red curves 
are exact solutions, black - approximations. 

by \j/{ri^) makes analysis easier. Another useful feature of (6.4) is the possibility to adjust the size of 
the low -frequency neighborhood of interest by changing Tj . 

To investigate quality of our approximations in the Fourier space, we analyze Fourier spectra of 
the exact Jacobian, velocity, stresses and their approximations. We consider only the second test case 
because the initial velocity in this case has full spectrum unlike the first test case. Since the approxima- 
tion tends to become worse at later times, we show the relevant spectra at the "worst case scenario" time 
t = 0.9. Fig. 1 1 depicts the spectra of the exact Jacobian and velocity (red curves) and their approxima- 
tions (black curves). The left panel of Fig. 1 1 indicates that we only capture about 70 first low frequency 
modes of the Jacobian. Similarly, the first 70 modes of the velocity are well reconstructed, while modes 
between 70 and 180 have much smaller amplitudes than in the exact velocity. While spectra of both 
Jacobian and velocity are not very accurate, spectral approximations of both convective and interaction 
stresses are quite good (see Fig. 12). Both exact stresses have only low frequency components (110 for 
the convective stress and only 70 for the interaction) and all these modes are captured perfectly! This 
demonstrates that it is not necessary to recover higher frequency modes of the Jacobian and velocity in 
order to approximate stress accurately. 




fcO.9 
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Fig. 12. Discrete Fourier coefficients of the convective stress (left panel) and interaction stress (right panel) for the second test 
case with N = \0, 000, B = 500, 77 = 0. 1 and t = 0.9. The red curves are exact solutions, black - approximations. 

Loss of accuracy in using truncated spectra for deconvolution often leads to Gibbs phenomenon. It 
is indeed present in both Jacobian and velocity reconstructions shown in Fig. 13. The amplitude of 
Gibbs ripples seem to increase with rj. In contrast, the approximations to stresses do not suffer from 
Gibb's oscillations as can be seen from Fig. 14. Gibb's phenomenon is typical for solutions of linear 
systems using a truncated SVD approach (see Boyd (2002); Bruno (2003); Bruno et al. (2007); Boyd 
& Ong (2009)). In Lyon (2012), Gibb's phenomenon is controlled by using Sobolev smoothing. In our 
case smoothing is done naturally by averaging present in the stress approximation (see (1.1)). 
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Fig. 14. Reconstruction of the convective stress (left panel) and interaction stress (right panel) in the second test case with 
N = 10,000, B = 500 and T] = 0.1. Exact and approximate solutions are shown at ^ = 0.9. The red curves are exact solutions, 
black - approximations. 



7. Scale separation with fixed t], B and varying 

In this section we investigate how the scale separation, i.e. ratio B to N affects the accuracy of recon- 
struction of the Jacobian and velocity. We use two test initial conditions as before with Gaussian window 
function 1/^(6) (x), 77 = 0. 1, B = 500 and = 1000, 2000, 5000 and 10,000. 

In the first test case, results of which are presented in Figs. 15 and 16, we observe that the error in 
reconstruction of the Jacobian (shown in the left panel of Fig. 15) increases as increases, while the 
error in velocity reconstruction (right panel of Fig. 15) does not depend on A^. Both errors oscillate in 
time and their oscillatory dynamics is related to the total computed energy oscillations depicted in Fig. 
7 in the left panel, i.e. when energy starts deviating from its initial (exact) value, error in the Jacobian 
reconstruction starts to increase and when the energy starts returning to its initial value, the error in the 
Jacobian starts decreasing. The error in the velocity approximation does not essentially depend on A^, 
oscillates as well and does not exceed 2% during the simulation time. Even though the reconstruction of 
the Jacobian becomes worse as the scale separation increases, the approximation of the convective stress 
(shown in the left panel of Fig. 16) gets better as A^ increases. At all times, the error in the convective 
stress does not exceed 3%. The error in the interaction stress does not depend on A^ as can be seen from 
the right panel of Fig. 16. 

In the second test case reconstruction of both Jacobian and velocity depends on A^ in a non- 
monotonic manner The intervals where the error with A^ = 10,000 is smaller and larger alternate in 
phase with oscillations of the total computed energy. While the error in the Jacobian approximation is 
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at most 0.2%, the error in the velocity approximation is rather large, especially at later times: it is below 
5% until t — 0.7 and then increases to 15 — 25% (15% for 10,000). The error in the approximation of 
the convective stress is shown in the left panel of Fig. 16. It is not monotonic in A^, and it does not 
exceed 6%. The error in the interaction stress does not depend on and stays below 2.5%. 
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Fig. 16. Effect of the scale separation on the error in approximation of the convective stress (left panel) and interaction stress 
(right panel) for the first test case with?] =0.1, v**' (,t), B = 500 and /V = 1000, 2000, 5000 and 10,000. 



8. Error estimates 

8.1 Estimates for filtered regularization methods 

In practice, the right hand side of (3.1) is known imprecisely, so instead of the exact b, one has an 
approximate vector b . The computed regularized solution is thus 

X =Lbj——^j. (8.1) 

.7=1 J 

Our goal is to estimate the error || x—x"'^ \\p, in some vector /7-norm. Usually p = 2, but here we 
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Fig. 17. Effect of the scale separation on the reconstruction of the Jacobian (left panel) and velocity (right panel) in the second 
test case with JJ =0.1, v/*"' (.r), B = 500 and W = 1000, 2000, 5000 and 10,000. 
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Fig. 18. Effect of the scale separation on the error in approximation of the convective stress (left panel) and interaction stress 
(right panel) for the second test case with 77 = 0. 1, v'*' [x), B = 500 and N = 1000, 2000, 5000 and 10,000. 



assume p E [1,°°)- By triangle inequality, 



X X 



a.S 



l-«/>(g,-,a) 

j=l 



D 

Lib J 



(8.2) 



v/=l 



The constant C{^,p) depends only on p and the components of the singular vectors ^j. 
Writing \bj\/Gj = \xj\, we see that the first term on the very right of (8.2) is bounded by 

Cilp)max\xj\ |l-0((j,,a)|'^y ^Ciitp) \\x\\^ ( f''^\l-(l>{f{t),a)\"dt]'' , (8.3) 







where we introduced a function /(f) : [Q,°°) — > (0, 1] that interpolates between singular values: 

/(;■) = (T;, /(0)-l, lm/(0-0. 

The function / is chosen to be continuous, non-negative, and strictly decreasing. To see that the integral 
is larger than the corresponding sum, note that the sum is the left-endpoint Riemann sum for the integral, 
and the function under the integral is increasing. 
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Similarly, 

Combining (8.2)-(8.4) we have 

\\x-x"^^\\p < Ci(|,;,)||x|U|l l-0(/(O,a)||L/'(o.D+i) (8.5) 
+ Ci(|,p) II b-b^ U\ <l>{m,OC)f-\t) ||i^(o,D+l) . 

By definition of (p, the first term can be made arbitrarily small by choosing a small enough. In 
the second term, as a — >^ 0, the norm of <p{f{t),a)f^^(t) typically increases. To control the second 
term, we need b — b to be small. This is typical of the error estimates available in the literature. Our 
inequalities differ from the standard ones because we use a /5-norm for the error, and oo-norms for x and 
b — b^. The standard estimates use 2-norms of x — x"'^ , b — b^, and what is essentially the oo-nomi 
for the (j)- and /-dependent terms. Depending on the actual (j) and /(f), our approach can yield tighter 
bounds. Improvement occurs if, loosely speaking, the integrals involving 0,/ in (8.3), (8.4) are smaller 
than maximal values of the integrands. 

Finally, we note that using Holder inequality with exponents q,q',q^^ + {q')^^ = 1 in the right hand 
side of (8.2) results in the estimates 

\\x-x"'^\\p < Czd, p,^) II a: llp^ll l-(^)(/(0, a) ||^,V(o,D+i) (8.6) 

+ C2{lp,q) II b-b' \\p,\\ mt),a)r\t) \\^,'^o,D+i) ■ 

8.2 Error in the interaction stress approximation 

The purpose of this section is to estimate the difference between the exact integral representation of the 
interaction stress T'l'.^^j in (2.14), and its closed form approximation 

Tln,){t,x) = ^^j V,ix-R) yu'{\p\)^Q„m{t,R+^-p)Qr,mit,R-lp)dp^ dR. 

(8.7) 

Since estimates will be local in time, we wiU suppress the dependence on t in the remainder of this 
section. Define the error 

m=Ti^){x)-Tint){^)- 

Next, introduce the abbreviated notation 

y+=y(/j+|p), j-=j{R-^-pi 

^{p) = U'{\p\)^. (8.8) 



and 

and denote 
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This function is smooth and can be assumed compactly supported on a shell D = {p : ci ^ |p| ^ C2} 
where ci > 0. With these notations, using an elementary identity 

aia2-bib2 =01(02-^2) +«2(«i -^1) - («i -^i)(«2-^2) 

we have 

7+7--e„[p'']+e^[p'']- = J+iJ--Qr^m~) + {J^-Qv\p'']^)QrilpT (8-9) 

= j+{j--Q,[p^]-)+riJ+-Q^m+) 

Now 

\E{x)\ < \Q\-hvip\xi/^\l l\0{p)\\J+J--Qr,m^Qn\p'^]~\{R,P)dpdR 
< |i2|~^sup|v^,,|sup|0| 

D 

J y'(^|/+ii7--e^[p'']-|+|7-||7+-e„[p'']+l+l7+-e^[p'']+ll7--e„[p'']-|)jpj/?. 

Changing variables in the last integral to 

and observing that the Jacobian of this transformation is and that the quantities marked by (respec- 
tively by ~) depend only on y-^ (respectively on ^2). we find 



2/ \J{yi)\dyJ \J-Qnm{y2)dy2+( f \J-Qnmiy)dy 



\E{x)\^e-''\Q\-^C{xjfr^,<P) 

Suppose now that / and Qrj [p^] are given by their discretizations on the fine mesh. Thus we can assume 
that they are piecewise constant functions having values Jj, Q-qlP^]] on the sets Sj C Q,j — 1,2,. ..,N 
of measure |^2|/A^. In this way, J and Gglp''] can be identified with, respectively, the vectors / = 

(71,72, • • • ,JnV and e = (Gi, 22, • • • , QnV- 

Suppose that there exists a constant M such that 

J^M. (8.10) 

With this, 



\E{x)\ < e-''|X2|-2c(VA^,0)(^yj (2M||J-e||i + ||J-e||?) (S.ll) 
= C(VA„,^)(2M||J-Q||i + ||J-Q||i). 

The last equality holds since N = £^ The norms are vector 1-norms that can be estimated using (8.6) 
with x = J, x"'^ = Q, and b representing a discretization of . 
The results of this section can be summarized in the 
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Theorem 8.1 Suppose that 

(i) J{t,x) satisfies (8.10) uniformly in /; 

(ii) defined in (8.8) is bounded; 

(iii) J{t,x) = Ji{t,x)Xj{x), Qri [p''] = LjLi Qi{t,x)Xj{x), where Xj are characteristic functions of 
sets 5^- such that Uj^i = Q., Sjr)Sk = <d if j^k, and \Sj\ =N-^\Q\. Then the error 



satisfies 



{int) {int) 



At) I < sup I I sup 1 0| (2M II J - e II 1 + II J - e II?) , 



where J = {Ji,J2, ■ ■ ■ jNf, G = (Gi, 62, • • • , QnV- 

The estimates for the error in T^^^ can be derived similarly, but would require more technical work 
because of the triple product structure of the integrand. It is also worth noting that, while a pointwise 
bound on J can be reasonably expected, similar bounds on the velocity v would blow up as e — )• 0. 
Consequently, estimating of the error T^^-^ — T^^^ is left to future work. 



9. Conclusions 

We study the numerical performance of the regularized deconvolution closiu'e introduced in Panchenko 
et al. (2011, subm). The closure method consists of the following. The average density and linear 
momentum are written as convolutions acting on respective fine scale functions: J and Jv, where J is 
the Jacobian of the inverse deformation map, and v is a particle velocity interpolant. These functions 
are approximately recovered by applying a regularized deconvolution to the averages. To construct 
the deconvolution operator, we use the theory of ill-posed problems. Closure is obtained by using 
these deconvolution approximations in the exact flux equations. This gives constitutive equations that 
express stress in terms of the average density and velocity. The exact stress is thus approximated by a 
sum of terms that have the "convolution sandwich" structure: they combine the convolution operator, 
a nonlinear composition or a product type operator, and the deconvolution operator. The resulting 
constitutive equations are nonlinear and nonlocal. 

The approximation quality depends on a choice of the window function y/ used to define averages, 
magnitude of scale separation, and values of the resolution and regularization parameters. Because of 
the nonlinearity of the problem, the error estimates tend to be too pessimistic. Therefore, we conduct 
numerical experiments to determine the dependence of the error on the above parameters. Since the 
Fourier spectrum of velocity seems to have a strong effect on the error, we consider two sets of initial 
conditions. In the first test case, the initial velocity is a low frequency mode sine function, while in the 
second test the initial velocity has full Fourier spectrum. The initial positions in both cases are equally 
spaced. 

We study window functions of different smoothness, varying from piecewise continuous to infinitely 
smooth. Among these functions, the Gaussian provides the best overall performance despite the fact that 
the corresponding integral deconvolution problem has the highest degree of ill-posedness. Numerical 
deconvolution amounts to solving an ill-conditioned linear system. We use a truncated SVD method 
with an additional spectral filtering of the right hand side. Filtering helps to reduce the effect of error 
that is present in every standard numerical SVD routine. 

The choice of the resolution parameter rj affects the size of the averaging window and the amount 
of high frequency filtering in the computed averages. Larger values of tj produce smoother and smaller 
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averages, thus causing the reconstruction to deteriorate. This tendency is counteracted by the presence of 
the convolution operator in the stress equations. We find that the overall error in the stress approximation 
tends to decrease with increasing 77 . Therefore, it is not necessary to have very good reconstruction of 
the Jacobian and velocity to have good approximations of the stresses. This "self-correcting" property 
is a noteworthy feature of the deconvolution closure. Another method to increase scale separation is to 
vary the number of particles while keeping rj fixed. The results in this case are less clear-cut compared 
to the case of increasing 77 . However, at times when the computed total energy is close to its exact value, 
the error in the stress decreases with increasing scale separation. 

The deconvolution error estimates derived in the paper are applicable to general SVD-based filtered 
regularization methods (see e.g. Kirsch (1996)). We also obtain error estimates for the interaction stress 
(the part of the total stress induced by interparticle forces). We believe that similar estimates can be also 
obtained for the remaining convective stress, but such estimates will be developed elsewhere. 

A. Window functions 

A window function xj/ is chosen to define a mesoscale average. This function has to satisfy several 
conditions: be nonnegative, fast decreasing, compactly supported (we also consider non-compactly 
supported functions like Gaussian), continuous and differentiable almost everywhere in the interior of 
its domain and /~ \j/{x) = 1. We use functions \jf^'\x), i = 1, ... ,6 of different order of smoothness 
starting from the characteristic function i/a(')(x) that is discontinuous atx = ±j up to infinitely many 
times differentiable Gaussian The window functions are defined in (A.1)-(A.6) and plotted in 



Fig. 2. 




(A.1) 




1/(2L), if \x\^L/2, 
-2/L{x-3L/2), if 
-2/L{x + 3L/2), if 
0, otherwise; 



L/2 < jc < 3L/2, 
-3L/2^x<-L/2, 



(A.2) 



-A/L}{x-L/2) 
A/L^{x+L/2), 



if 



0<x<L/2, 
-L/2<jc<0, 



Va(3)(jc) = < 



if 



(A.3) 



0, otherwise; 




(A.4) 




(A.5) 




(A.6) 
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B. Lennard-Jones potential 

The dynamics of particles considered in this paper is governed by Lennard-Jones potential defined as 



with the potential well depth e = 0.025. The same potential but with e = 0.25 was used in Panchenko 
et al. (subm). The magnitude of £ defines how strong interaction between particles is. The potential 
is zero at the distance given by a and reaches its minimum at the distance h = 2^/^(7 at which parti- 
cles are in equilibrium. For smaller distances ^ < h, the potential is repulsive whereas for ^ > it is 
mildly attractive. When the distance § > 2.5h, the force is very small and we set it to zero to speed up 
computations. This truncation of the potential tail typically takes into account 3 particles on each side 
from a current particle. Truncating at larger distances slightly decreases deviations of the total computed 
energy from its exact value (at t = 0) and as a result slightly decreases the error in the approximation of 
the stresses, more so in the interaction stress. 
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